rm(list = ls(all = TRUE))
gc()
##           used (Mb) gc trigger (Mb) max used (Mb)
## Ncells  555542 29.7    1244034 66.5   686457 36.7
## Vcells 1022955  7.9    8388608 64.0  1876634 14.4

1 packages

# %>% 
library(magrittr)
library(ggplot2)


`%!in%` = Negate(`%in%`)

2 palette

n = 4
col = RColorBrewer::brewer.pal(n, 'Paired')

pal = col[seq(1, n, 2)]
my_dir = file.path("..", "reports", "Table7")
if (!dir.exists(my_dir)) {
  dir.create(my_dir)
}
fr = file.path(my_dir, 'Tests.txt')
file.create(fr)
## [1] TRUE

3 f-ctions

3.1 summary stat

summary_stats <- function(df, measurevar, groupvars, conf.level = 0.95) {
  df %>%
    dplyr::group_by(across(all_of(groupvars))) %>%
    dplyr::summarise(
      N = sum(!is.na(.data[[measurevar]])),
      mean = mean(.data[[measurevar]], na.rm = TRUE),
      sd = sd(.data[[measurevar]], na.rm = TRUE),
      .groups = "drop"
    ) %>%
    dplyr::mutate(
      se = sd / sqrt(N),
      ci = se * qt(conf.level/2 + 0.5, N - 1)
    )
}

3.2 scale PS

process_data <- function(df, groupvars, measurevar, scale_treatment = "noninoculated") {

  df.long = df %>%
    tidyr::pivot_longer(cols = 3:ncol(df), names_to = "transcript",
                        values_to = measurevar, values_drop_na = TRUE)
  data.SE = summary_stats(df.long, measurevar, groupvars)
  scale_reference = data.SE %>%
    dplyr::filter(Treatment == scale_treatment) %>%
    dplyr::select(Tissue, transcript, mean) %>%
    dplyr::rename(scale_mean = mean)
  df.scaled <- dplyr::left_join(df.long, scale_reference, by = c("Tissue", "transcript")) %>%
    dplyr::mutate(scaled = .data[[measurevar]] / scale_mean) %>%
    dplyr::select(-scale_mean)
  df.scaled = df.scaled %>% 
    dplyr::arrange(dplyr::desc(transcript), dplyr::desc(Treatment))
  shoots = df.scaled %>% dplyr::filter(Tissue == "shoots")
  roots = df.scaled %>% dplyr::filter(Tissue == "roots")
  
  return(list(
    shoots = shoots,
    roots = roots
  ))
}

3.3 perm t

perm_test_by_transcript <- function(mydata.long, measurevar = "measurement", groupvar = "Treatment") {
  
  temp = data.frame()
  transcript_levels = levels(mydata.long$transcript)
  treatment_levels = levels(mydata.long[[groupvar]])
  treatment_pairs = combn(treatment_levels, 2, simplify = FALSE)
  
  for(i in transcript_levels) {
    data = mydata.long[mydata.long$transcript == i, ]
    k = 8 # nrow(data)
    pvalues = purrr::map_dbl(treatment_pairs, function(x) {
      subset_data = base::subset(data, data[[groupvar]] %in% x)
      res = MKinfer::perm.t.test(
        formula = stats::as.formula(base::paste(measurevar, "~", groupvar)),
        data = subset_data,
        alternative = "two.sided",
        mu = 0,
        paired = FALSE,
        var.equal = FALSE,
        conf.level = 0.95,
        perm.conf.int = 0.95,
        symmetric = TRUE,
        p.adjust.method = "holm",
        detailed = TRUE,
        R = choose(k, k/2), # sum(choose(k, 1:(k-1))),
        set.seed = 123456
      )
      res$perm.p.value
    })
    
    tmp = base::as.data.frame(base::t(pvalues))
    colnames(tmp) = purrr::map_chr(treatment_pairs, ~base::paste(.x, collapse = ' vs '))
    rownames(tmp) = i
    
    temp <- base::rbind(temp, tmp)
  }
  
  return(temp)
}

3.4 plot

plot_gene_expression <- function(data.SE
                                 , data.long
                                 , stat.test.sig
                                 , transcripts_excl = c('13-LOX', 'PTI5')
                                 , color_values
                                 , facet_cols = 6
                                 , y_scales = NULL
                                 , dodge_width = 0.8,
                                 plot_title = ""
                                 ) {
  
  dodge = position_dodge(width = dodge_width)
  
  data.SE.filtered = dplyr::filter(data.SE, !transcript %in% transcripts_excl)
  data.long.filtered = dplyr::filter(data.long, !transcript %in% transcripts_excl)
  
  stat.test.filtered = NULL
  if (nrow(stat.test.sig) > 0) stat.test.filtered = dplyr::filter(stat.test.sig, !transcript %in% transcripts_excl)
  
  p = ggplot(data.SE.filtered, aes(x = Treatment, y = mean), colour = "black") +
    geom_point(size=3.5, shape = 22, position = dodge, aes(fill = Treatment), colour = "black") +
    geom_point(data = data.long.filtered, aes(x = Treatment, y = scaled, fill = Treatment), colour = "black",
               size = 2.0, shape = 21,
               position = dodge) +
    geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.3, lwd = 0.5,
                  position = dodge) +
    facet_wrap(~ transcript, ncol = facet_cols, scales = "free", drop = TRUE)
  
  if (!is.null(y_scales)) {
    p = p + ggh4x::facetted_pos_scales(y = y_scales)
  }
  
  p = p +
    # geom_hline(yintercept = 2.5, alpha = 0.0) +
    geom_hline(yintercept = 0, alpha = 0.0) +
    geom_hline(yintercept = 1, alpha = 0.5, linetype = "dotted", col = "gray45") +
    ggtitle(plot_title) +
    theme_bw() +
    scale_colour_manual(name = "", values = rev(color_values)) +
    scale_fill_manual(name = "", values = rev(color_values)) +
    labs(x = "", y = "Relative gene expression (+/- SE)"
         # , subtitle = "Permutation t-test measurements"
         ) +
    theme(plot.subtitle = element_text(size = 10),
          axis.text = element_text(size = 12.5),
          axis.text.x = element_text(size = 12.5, angle = 90),
          axis.title = element_text(size = 12.5, face = "bold"),
          strip.text = element_text(size = 12.5),
          title = element_text(size = 15, face = "bold"),
          # axis.ticks.x = element_blank(),
          legend.key.height = unit(1.5, "cm"),
          legend.key.width = unit(1.75, "cm"),
          legend.text = element_text(size = 12.5),
          legend.title = element_text(size = 12.5),
          legend.background = element_rect(fill = "transparent", size = 0.5, linetype = "dotted"),
          legend.position = "top",
          legend.justification = "right",
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.title.y.right = element_blank(),
          axis.text.y.right = element_blank(),
          # axis.ticks.y = element_blank(),
          axis.text.y = element_text(size = 12.5, margin = margin(r = 0)),
          panel.spacing = unit(1, "lines"),
          strip.background = element_rect(size = 0.5, fill = "transparent", color = NA) ,
          panel.border = element_blank(),
          axis.line = element_line(color = "black")
          )
  
    p = p + theme(legend.position = "none")

  


  
  if (!is.null(stat.test.filtered) && nrow(stat.test.filtered) > 0) {
    
    p = p + ggpubr::stat_pvalue_manual(
      stat.test.filtered,
      label = "p.adj.signif",
      xmin = "xmin",
      xmax = "xmax",
      y.position = "y.position",
      hide.ns = FALSE,
      tip.length = 0.01,
      step.increase = 0.075,
      inherit.aes = FALSE 
    )
    
  }
    
  return(p)
  
}

3.5 effects

Denote: PR1B has LOQ values

dens_and_effect <- function(mydata.long, p_colors) {
  
  # Plot density with x="measurement"
  p1 = ggpubr::ggdensity(mydata.long, x = "measurement",
              add = "mean", rug = TRUE,
              color = "Treatment", fill = "Treatment",
              facet.by = 'transcript') +
    scale_fill_manual(name = "", values = rev(p_colors)) +
    scale_color_manual(name = "", values = rev(p_colors)) +
    ggtitle("measurement") +
    facet_wrap(~transcript, ncol = 4, scales = "free")
  print(p1)
  
  # Plot density with x="scaled"
  p2 = ggpubr::ggdensity(mydata.long, x = "scaled",
              add = "mean", rug = TRUE,
              color = "Treatment", fill = "Treatment",
              facet.by = 'transcript') +
    scale_fill_manual(name = "", values = rev(p_colors)) +
    scale_color_manual(name = "", values = rev(p_colors)) +
    ggtitle("scaled") +
    facet_wrap(~transcript, ncol = 4, scales = "free")
  print(p2)
  
  
  cat(crayon::red('Shapiro-Wilk Normality Test\n'))
  
  # Shapiro-Wilk normality test on measurement
  norm_meas = mydata.long %>%
    dplyr::select(-Tissue, -Treatment, -scaled) %>%
    tidyr::gather(key = "transcript", value = "measurement") %>%
    dplyr::group_by(transcript) %>%
    dplyr::do(broom::tidy(shapiro.test(.$measurement))) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(p.adj = p.adjust(p.value, method = "holm")) %>%  # or "bonferroni", "BH"
    dplyr::select(transcript, W = statistic, p.value, p.adj)
  print(norm_meas)
  
  # Shapiro-Wilk normality test on scaled
  norm_scaled = mydata.long %>%
    dplyr::select(-Tissue, -Treatment, -measurement) %>%
    tidyr::gather(key = "transcript", value = "scaled") %>%
    dplyr::group_by(transcript) %>%
    dplyr::do(broom::tidy(shapiro.test(.$scaled))) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(p.adj = p.adjust(p.value, method = "holm")) %>%  # or "bonferroni", "BH"
    dplyr::select(transcript, W = statistic, p.value, p.adj)
  print(norm_scaled)
  

  cat(crayon::red("Quantile-Quantile plots\n"))
  plot_list = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      ggpubr::ggqqplot(.x$measurement) +
        ggplot2::ggtitle(paste("", .y$transcript))
    })
  ggpubr::ggarrange(plotlist = plot_list, ncol = 4, nrow = ceiling(length(plot_list) / 4))

  
  cat(crayon::red("Levene's test for homogeneity of variance across groups\n"))
  
  # Levene's test on measurement
  lev_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::levene_test(measurement ~ Treatment)
  print(lev_meas)
  
  # Levene's test on scaled
  lev_scaled = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::levene_test(scaled ~ Treatment)
  print(lev_scaled)
  
  cat(crayon::red("Wilcoxon effect size\n"))
  
  # Wilcoxon effect size on measurement
  eff_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::wilcox_effsize(measurement ~ Treatment)
  print(eff_meas)
  
  # Wilcoxon effect size on scaled
  eff_scaled = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::wilcox_effsize(scaled ~ Treatment)
  print(eff_scaled)
  
  cat(crayon::red("Cohen's d Measure of Effect Size\n"))
  
  # Cohen's d on measurement
  coh_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::cohens_d(measurement ~ Treatment, paired = FALSE)
  print(coh_meas)
  
  # Cohen's d on scaled
  coh_scaled = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::cohens_d(scaled ~ Treatment, paired = FALSE)
  print(coh_scaled)
  
  invisible(list(
    p_scaled = p1,
    p_measurement = p2,
    shapiro_measurement = norm_meas,
    shapiro_scaled = norm_scaled,
    levene_measurement = lev_meas,
    levene_scaled = lev_scaled,
    wilcoxon_measurement = eff_meas,
    wilcoxon_scaled = eff_scaled,
    cohensd_measurement = coh_meas,
    cohensd_scaled = coh_scaled
  ))
}

3.6 test and plot

test_and_plot <- function(data_long_raw
                          , myorder
                          , pal
                          , what
                          , plot_gene_expression_func
                          , groupvars = c("Treatment", "transcript")
                          , y_scales1
                          , y_scales2
                          ) {

  
  mydata.long = dplyr::as_tibble(data.table::data.table(data_long_raw))
  mydata.long$transcript = factor(mydata.long$transcript, levels = myorder)
  mydata.long = mydata.long %>% dplyr::arrange(factor(transcript, levels = myorder))
  
  mydata.long.SE = summary_stats(mydata.long, measurevar = "scaled", groupvars = groupvars)
  
  results = dens_and_effect(mydata.long, p_colors = pal)
  
  cat(what, file = fr, append = TRUE, sep = "\n")
  cat("shapiro", file = fr, append = TRUE, sep = "\n")
  output_text = capture.output(as.data.frame(results$shapiro_measurement))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("levene", file = fr, append = TRUE, sep = "\n")
  output_text = capture.output(as.data.frame(results$levene_measurement))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("wilcoxon", file = fr, append = TRUE, sep = "\n")
  output_text = capture.output(as.data.frame(results$wilcoxon_measurement))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("cohensd", file = fr, append = TRUE, sep = "\n")
  output_text = capture.output(as.data.frame(results$cohensd_measurement))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")

  
  temp = perm_test_by_transcript(mydata.long, measurevar = "measurement", groupvar = "Treatment")
  temp$transcript = rownames(temp)
  perm = tidyr::gather(temp, contrast, perm.p, colnames(temp)[1]:colnames(temp)[ncol(temp)-1], factor_key=TRUE)
  perm$group1 = gsub(' vs.*', '', perm$contrast)
  perm$group2 = gsub('.* vs ', '', perm$contrast)
  perm$perm.p.adj = p.adjust(perm$perm.p, method = 'BH')
  perm$perm.p.adj.signif = 'ns'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.0001] = '**'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.001] = '**'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.05] = '*'
  
  # like rstatix object
  stat.test = perm %>%
    dplyr::select(transcript, group1, group2, perm.p, perm.p.adj, perm.p.adj.signif) %>%
    dplyr::rename(
      p = perm.p,
      p.adj = perm.p.adj,
      p.adj.signif = perm.p.adj.signif
    )
  y_pos_df = mydata.long %>%
    dplyr::group_by(transcript, Treatment) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE)) %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(y.position = max(max_scaled) * 1.0)
  stat.test = stat.test %>%
    dplyr::left_join(y_pos_df, by = "transcript")
  stat.test = stat.test %>%
    dplyr::mutate(
      xmin = as.numeric(factor(group1, levels = unique(mydata.long$Treatment))),
      xmax = as.numeric(factor(group2, levels = unique(mydata.long$Treatment)))
    )
  stat.test.sig <- stat.test %>%
    dplyr::filter(p.adj < 0.05)
  
  group1 =  sapply(y_scales1, function(x) {
    s = deparse(x)
    s_full = paste(s, collapse = " ")
    sub('.*transcript == *"([^"]+)".*', '\\1', s_full)
  })
  group2 =  sapply(y_scales2, function(x) {
    s = deparse(x)
    s_full = paste(s, collapse = " ")
    sub('.*transcript == *"([^"]+)".*', '\\1', s_full)
  })


  p1 = plot_gene_expression_func(
    data.SE = mydata.long.SE,
    data.long = mydata.long,
    stat.test.sig = stat.test.sig,
    transcripts_excl = group2,
    facet_cols = 2,
    color_values = pal,
    plot_title = what,
    y_scales = y_scales1,
    dodge_width = 0.8
  )
  print(p1)
  
  p2 = plot_gene_expression_func(
    data.SE = mydata.long.SE,
    data.long = mydata.long,
    stat.test.sig = stat.test.sig,
    transcripts_excl = group1,
    facet_cols = 6,
    color_values = pal,
    plot_title = what,
    y_scales = y_scales2,
    dodge_width = 0.8
  )
  print(p2)
  
  return(list(plot1 = p1, plot2 = p2, stat.test = stat.test))
}

4 Test

4.1 input

fp = file.path('..', "input")
list.files(fp)
## [1] "README.md"     "Table S1.xlsx" "Table S2.xlsx" "Table S7.xlsx"
fn = 'Table S7.xlsx'

PS216 = openxlsx::read.xlsx(xlsxFile = file.path(fp, fn),
                            sheet = 'Test',
                            startRow = 1,
                            colNames = TRUE,
                            rowNames = FALSE,
                            detectDates = FALSE,
                            skipEmptyRows = TRUE,
                            skipEmptyCols = TRUE,
                            rows = NULL,
                            cols = NULL,
                            check.names = FALSE,
                            sep.names = ".",
                            namedRegion = NULL,
                            na.strings = "NA",
                            fillMergedCells = FALSE)

data.table::setDT(PS216)
PS216[, Sample.ID := NULL]
PS216[, Genotype := NULL]
PS216$Tissue = as.factor(trimws(PS216$Tissue))
PS216$Treatment = factor(trimws(PS216$Treatment), levels = c("noninoculated", "inoculated"))

PS218 = PS216[grep('PS-216', PS216$Strain, invert = TRUE), ]
PS216 = PS216[grep('PS-218', PS216$Strain, invert = TRUE), ]
PS216[, Strain := NULL]
PS218[, Strain := NULL]

4.2 params

groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("RBOHD", "PR1B", "CPI8", "CAB", "BGLU2", "HSP70", "PTI5", "13-LOX")

4.3 results

shoots.216 = process_data(PS216, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
roots.216 = process_data(PS216, groupvars, measurevar, scale_treatment = "noninoculated")$roots
shoots.218 = process_data(PS218, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
roots.218 = process_data(PS218, groupvars, measurevar, scale_treatment = "noninoculated")$roots


y_scales1 = list(transcript == "PTI5" ~ scale_y_continuous(limits = c(0, 11), breaks = seq(0, 11, 1))
                 , transcript == "13-LOX" ~ scale_y_continuous(limits = c(0, 11), breaks = seq(0, 11, 1)))
y_scales2 = list(transcript == "RBOHD" ~ scale_y_continuous(limits = c(0, NA), breaks = seq(0, 11, 0.5))
                 , transcript == "PR1B" ~ scale_y_continuous(limits = c(0,NA), breaks = seq(0, 11, 1))
                 , transcript == "CPI8" ~ scale_y_continuous(limits = c(0, NA), breaks = seq(0, 11, 1))
                 , transcript == "CAB" ~ scale_y_continuous(limits = c(0, NA), breaks = seq(0, 11, 1))
                 , transcript == "BGLU2" ~ scale_y_continuous(limits = c(0, NA), breaks = seq(0, 11, 0.5))
                 , transcript == "HSP70" ~ scale_y_continuous(limits = c(0, NA), breaks = seq(0, 11, 0.5))
                 )
result = test_and_plot(data_long_raw = shoots.216
                       , myorder = myorder
                       , pal = pal
                       , what = "shoots 216"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## Shapiro-Wilk Normality Test
## # A tibble: 8 × 4
##   transcript     W  p.value   p.adj
##   <fct>      <dbl>    <dbl>   <dbl>
## 1 RBOHD      0.947 0.682    1      
## 2 PR1B       0.601 0.000161 0.00129
## 3 CPI8       0.918 0.415    1      
## 4 CAB        0.953 0.745    1      
## 5 BGLU2      0.930 0.520    1      
## 6 HSP70      0.962 0.829    1      
## 7 PTI5       0.879 0.186    1      
## 8 13-LOX     0.790 0.0226   0.158  
## # A tibble: 8 × 4
##   transcript     W  p.value   p.adj
##   <fct>      <dbl>    <dbl>   <dbl>
## 1 RBOHD      0.947 0.682    1      
## 2 PR1B       0.601 0.000161 0.00129
## 3 CPI8       0.918 0.415    1      
## 4 CAB        0.953 0.745    1      
## 5 BGLU2      0.930 0.520    1      
## 6 HSP70      0.962 0.829    1      
## 7 PTI5       0.879 0.186    1      
## 8 13-LOX     0.790 0.0226   0.158  
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6   0.181   0.685
## 2 PR1B           1     6   0.478   0.515
## 3 CPI8           1     6   0.725   0.427
## 4 CAB            1     6   0.00927 0.926
## 5 BGLU2          1     6   0.739   0.423
## 6 HSP70          1     6   0.0138  0.910
## 7 PTI5           1     6   1.69    0.242
## 8 13-LOX         1     6   2.61    0.157
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6   0.181   0.685
## 2 PR1B           1     6   0.478   0.515
## 3 CPI8           1     6   0.725   0.427
## 4 CAB            1     6   0.00927 0.926
## 5 BGLU2          1     6   0.739   0.423
## 6 HSP70          1     6   0.0138  0.910
## 7 PTI5           1     6   1.69    0.242
## 8 13-LOX         1     6   2.61    0.157
## Wilcoxon effect size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.102 RBOHD          4     4 small    
## 2 measurement noninoculated inoculated   0.306 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated   0.510 CPI8           4     4 large    
## 4 measurement noninoculated inoculated   0.102 CAB            4     4 small    
## 5 measurement noninoculated inoculated   0.612 BGLU2          4     4 large    
## 6 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.714 13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.102 RBOHD          4     4 small    
## 2 scaled noninoculated inoculated   0.306 PR1B           4     4 moderate 
## 3 scaled noninoculated inoculated   0.510 CPI8           4     4 large    
## 4 scaled noninoculated inoculated   0.102 CAB            4     4 small    
## 5 scaled noninoculated inoculated   0.612 BGLU2          4     4 large    
## 6 scaled noninoculated inoculated   0.204 HSP70          4     4 small    
## 7 scaled noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 scaled noninoculated inoculated   0.714 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.269 RBOHD          4     4 small    
## 2 measurement noninoculated inoculated  -0.696 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated  -1.21  CPI8           4     4 large    
## 4 measurement noninoculated inoculated  -0.425 CAB            4     4 small    
## 5 measurement noninoculated inoculated  -1.44  BGLU2          4     4 large    
## 6 measurement noninoculated inoculated  -0.389 HSP70          4     4 small    
## 7 measurement noninoculated inoculated  -4.45  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -1.68  13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated  -0.269 RBOHD          4     4 small    
## 2 scaled noninoculated inoculated  -0.696 PR1B           4     4 moderate 
## 3 scaled noninoculated inoculated  -1.21  CPI8           4     4 large    
## 4 scaled noninoculated inoculated  -0.425 CAB            4     4 small    
## 5 scaled noninoculated inoculated  -1.44  BGLU2          4     4 large    
## 6 scaled noninoculated inoculated  -0.389 HSP70          4     4 small    
## 7 scaled noninoculated inoculated  -4.45  PTI5           4     4 large    
## 8 scaled noninoculated inoculated  -1.68  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p", colnames(result$stat.test))]
print(res)
##   transcript        group1     group2          p     p.adj p.adj.signif
## 1      RBOHD noninoculated inoculated 0.54285714 0.5857143           ns
## 2       PR1B noninoculated inoculated 0.40000000 0.5857143           ns
## 3       CPI8 noninoculated inoculated 0.18571429 0.3714286           ns
## 4        CAB noninoculated inoculated 0.58571429 0.5857143           ns
## 5      BGLU2 noninoculated inoculated 0.02857143 0.1142857           ns
## 6      HSP70 noninoculated inoculated 0.51428571 0.5857143           ns
## 7       PTI5 noninoculated inoculated 0.01428571 0.1142857           ns
## 8     13-LOX noninoculated inoculated 0.04285714 0.1142857           ns
output_text = capture.output(paste0("Test - shoots - 216"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
output_text = capture.output(res)
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")
  

ggsave(
  filename = file.path(my_dir, "Test_shoots.216_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_shoots.216_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)



result = test_and_plot(data_long_raw = roots.216
                       , myorder = myorder
                       , pal = pal
                       , what = "roots 216"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## Shapiro-Wilk Normality Test
## # A tibble: 8 × 4
##   transcript     W   p.value    p.adj
##   <fct>      <dbl>     <dbl>    <dbl>
## 1 RBOHD      0.921 0.438     1       
## 2 PR1B       0.528 0.0000222 0.000177
## 3 CPI8       0.824 0.0513    0.257   
## 4 CAB        0.786 0.0204    0.122   
## 5 BGLU2      0.555 0.0000471 0.000330
## 6 HSP70      0.955 0.765     1       
## 7 PTI5       0.828 0.0572    0.257   
## 8 13-LOX     0.915 0.391     1       
## # A tibble: 8 × 4
##   transcript     W   p.value    p.adj
##   <fct>      <dbl>     <dbl>    <dbl>
## 1 RBOHD      0.921 0.438     1       
## 2 PR1B       0.528 0.0000222 0.000177
## 3 CPI8       0.824 0.0513    0.257   
## 4 CAB        0.786 0.0204    0.122   
## 5 BGLU2      0.555 0.0000471 0.000330
## 6 HSP70      0.955 0.765     1       
## 7 PTI5       0.828 0.0572    0.257   
## 8 13-LOX     0.915 0.391     1       
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 RBOHD          1     6    1.19   0.317 
## 2 PR1B           1     6    1.00   0.356 
## 3 CPI8           1     6    7.39   0.0347
## 4 CAB            1     6    3.90   0.0956
## 5 BGLU2          1     6    0.792  0.408 
## 6 HSP70          1     6    0.0123 0.915 
## 7 PTI5           1     6    9.03   0.0238
## 8 13-LOX         1     6    0.258  0.630 
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 RBOHD          1     6    1.19   0.317 
## 2 PR1B           1     6    1.00   0.356 
## 3 CPI8           1     6    7.39   0.0347
## 4 CAB            1     6    3.90   0.0956
## 5 BGLU2          1     6    0.792  0.408 
## 6 HSP70          1     6    0.0123 0.915 
## 7 PTI5           1     6    9.03   0.0238
## 8 13-LOX         1     6    0.258  0.630 
## Wilcoxon effect size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.306 RBOHD          4     4 moderate 
## 2 measurement noninoculated inoculated   0.102 PR1B           4     4 small    
## 3 measurement noninoculated inoculated   0.408 CPI8           4     4 moderate 
## 4 measurement noninoculated inoculated   0.408 CAB            4     4 moderate 
## 5 measurement noninoculated inoculated   0.306 BGLU2          4     4 moderate 
## 6 measurement noninoculated inoculated   0.816 HSP70          4     4 large    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.306 RBOHD          4     4 moderate 
## 2 scaled noninoculated inoculated   0.102 PR1B           4     4 small    
## 3 scaled noninoculated inoculated   0.408 CPI8           4     4 moderate 
## 4 scaled noninoculated inoculated   0.408 CAB            4     4 moderate 
## 5 scaled noninoculated inoculated   0.306 BGLU2          4     4 moderate 
## 6 scaled noninoculated inoculated   0.816 HSP70          4     4 large    
## 7 scaled noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 scaled noninoculated inoculated   0.816 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.318 RBOHD          4     4 small    
## 2 measurement noninoculated inoculated  -0.653 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated  -1.26  CPI8           4     4 large    
## 4 measurement noninoculated inoculated  -1.60  CAB            4     4 large    
## 5 measurement noninoculated inoculated  -0.499 BGLU2          4     4 small    
## 6 measurement noninoculated inoculated  -2.82  HSP70          4     4 large    
## 7 measurement noninoculated inoculated  -2.43  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -3.09  13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.318 RBOHD          4     4 small    
## 2 scaled noninoculated inoculated  -0.653 PR1B           4     4 moderate 
## 3 scaled noninoculated inoculated  -1.26  CPI8           4     4 large    
## 4 scaled noninoculated inoculated  -1.60  CAB            4     4 large    
## 5 scaled noninoculated inoculated  -0.499 BGLU2          4     4 small    
## 6 scaled noninoculated inoculated  -2.82  HSP70          4     4 large    
## 7 scaled noninoculated inoculated  -2.43  PTI5           4     4 large    
## 8 scaled noninoculated inoculated  -3.09  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p", colnames(result$stat.test))]
print(res)
##   transcript        group1     group2          p      p.adj p.adj.signif
## 1      RBOHD noninoculated inoculated 0.54285714 0.62040816           ns
## 2       PR1B noninoculated inoculated 0.47142857 0.62040816           ns
## 3       CPI8 noninoculated inoculated 0.14285714 0.22857143           ns
## 4        CAB noninoculated inoculated 0.10000000 0.20000000           ns
## 5      BGLU2 noninoculated inoculated 0.94285714 0.94285714           ns
## 6      HSP70 noninoculated inoculated 0.01428571 0.03809524            *
## 7       PTI5 noninoculated inoculated 0.01428571 0.03809524            *
## 8     13-LOX noninoculated inoculated 0.01428571 0.03809524            *
output_text = capture.output(paste0("Test - roots - 216"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
output_text = capture.output(res)
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")


ggsave(
  filename = file.path(my_dir, "Test_roots.216_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_roots.216_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)



result = test_and_plot(data_long_raw = shoots.218
                       , myorder = myorder
                       , pal = pal
                       , what = "shoots 218"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## Shapiro-Wilk Normality Test
## # A tibble: 8 × 4
##   transcript     W   p.value    p.adj
##   <fct>      <dbl>     <dbl>    <dbl>
## 1 RBOHD      0.943 0.644     1       
## 2 PR1B       0.561 0.0000546 0.000437
## 3 CPI8       0.744 0.00700   0.0490  
## 4 CAB        0.945 0.658     1       
## 5 BGLU2      0.906 0.326     1       
## 6 HSP70      0.939 0.598     1       
## 7 PTI5       0.906 0.326     1       
## 8 13-LOX     0.898 0.279     1       
## # A tibble: 8 × 4
##   transcript     W   p.value    p.adj
##   <fct>      <dbl>     <dbl>    <dbl>
## 1 RBOHD      0.943 0.644     1       
## 2 PR1B       0.561 0.0000546 0.000437
## 3 CPI8       0.744 0.00700   0.0490  
## 4 CAB        0.945 0.658     1       
## 5 BGLU2      0.906 0.326     1       
## 6 HSP70      0.939 0.598     1       
## 7 PTI5       0.906 0.326     1       
## 8 13-LOX     0.898 0.279     1       
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6  0.255    0.631
## 2 PR1B           1     6  0.585    0.473
## 3 CPI8           1     6  3.74     0.101
## 4 CAB            1     6  0.0364   0.855
## 5 BGLU2          1     6  0.00507  0.946
## 6 HSP70          1     6  0.000737 0.979
## 7 PTI5           1     6  0.0321   0.864
## 8 13-LOX         1     6  1.91     0.216
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6  0.255    0.631
## 2 PR1B           1     6  0.585    0.473
## 3 CPI8           1     6  3.74     0.101
## 4 CAB            1     6  0.0364   0.855
## 5 BGLU2          1     6  0.00507  0.946
## 6 HSP70          1     6  0.000737 0.979
## 7 PTI5           1     6  0.0321   0.864
## 8 13-LOX         1     6  1.91     0.216
## Wilcoxon effect size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.408 RBOHD          4     4 moderate 
## 2 measurement noninoculated inoculated   0.204 PR1B           4     4 small    
## 3 measurement noninoculated inoculated   0     CPI8           4     4 small    
## 4 measurement noninoculated inoculated   0.306 CAB            4     4 moderate 
## 5 measurement noninoculated inoculated   0.612 BGLU2          4     4 large    
## 6 measurement noninoculated inoculated   0.612 HSP70          4     4 large    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.612 13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.408 RBOHD          4     4 moderate 
## 2 scaled noninoculated inoculated   0.204 PR1B           4     4 small    
## 3 scaled noninoculated inoculated   0     CPI8           4     4 small    
## 4 scaled noninoculated inoculated   0.306 CAB            4     4 moderate 
## 5 scaled noninoculated inoculated   0.612 BGLU2          4     4 large    
## 6 scaled noninoculated inoculated   0.612 HSP70          4     4 large    
## 7 scaled noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 scaled noninoculated inoculated   0.612 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -1.23  RBOHD          4     4 large    
## 2 measurement noninoculated inoculated  -0.702 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated   0.436 CPI8           4     4 small    
## 4 measurement noninoculated inoculated  -0.615 CAB            4     4 moderate 
## 5 measurement noninoculated inoculated  -1.79  BGLU2          4     4 large    
## 6 measurement noninoculated inoculated  -1.77  HSP70          4     4 large    
## 7 measurement noninoculated inoculated  -4.62  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -1.37  13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated  -1.23  RBOHD          4     4 large    
## 2 scaled noninoculated inoculated  -0.702 PR1B           4     4 moderate 
## 3 scaled noninoculated inoculated   0.436 CPI8           4     4 small    
## 4 scaled noninoculated inoculated  -0.615 CAB            4     4 moderate 
## 5 scaled noninoculated inoculated  -1.79  BGLU2          4     4 large    
## 6 scaled noninoculated inoculated  -1.77  HSP70          4     4 large    
## 7 scaled noninoculated inoculated  -4.62  PTI5           4     4 large    
## 8 scaled noninoculated inoculated  -1.37  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p", colnames(result$stat.test))]
print(res)
##   transcript        group1     group2          p     p.adj p.adj.signif
## 1      RBOHD noninoculated inoculated 0.11428571 0.1828571           ns
## 2       PR1B noninoculated inoculated 0.37142857 0.4952381           ns
## 3       CPI8 noninoculated inoculated 0.70000000 0.7000000           ns
## 4        CAB noninoculated inoculated 0.48571429 0.5551020           ns
## 5      BGLU2 noninoculated inoculated 0.05714286 0.1828571           ns
## 6      HSP70 noninoculated inoculated 0.11428571 0.1828571           ns
## 7       PTI5 noninoculated inoculated 0.08571429 0.1828571           ns
## 8     13-LOX noninoculated inoculated 0.05714286 0.1828571           ns
output_text = capture.output(paste0("Test - shoots - 218"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
output_text = capture.output(res)
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")



ggsave(
  filename = file.path(my_dir, "Test_shoots.218_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_shoots.218_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)


result = test_and_plot(data_long_raw = roots.218
                       , myorder = myorder
                       , pal = pal
                       , what = "roots 218"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## Shapiro-Wilk Normality Test
## # A tibble: 8 × 4
##   transcript     W  p.value   p.adj
##   <fct>      <dbl>    <dbl>   <dbl>
## 1 RBOHD      0.875 0.170    0.678  
## 2 PR1B       0.637 0.000424 0.00339
## 3 CPI8       0.802 0.0300   0.180  
## 4 CAB        0.777 0.0162   0.114  
## 5 BGLU2      0.894 0.254    0.763  
## 6 HSP70      0.899 0.280    0.763  
## 7 PTI5       0.826 0.0534   0.267  
## 8 13-LOX     0.906 0.327    0.763  
## # A tibble: 8 × 4
##   transcript     W  p.value   p.adj
##   <fct>      <dbl>    <dbl>   <dbl>
## 1 RBOHD      0.875 0.170    0.678  
## 2 PR1B       0.637 0.000424 0.00339
## 3 CPI8       0.802 0.0300   0.180  
## 4 CAB        0.777 0.0162   0.114  
## 5 BGLU2      0.894 0.254    0.763  
## 6 HSP70      0.899 0.280    0.763  
## 7 PTI5       0.826 0.0534   0.267  
## 8 13-LOX     0.906 0.327    0.763  
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 RBOHD          1     6    8.01   0.0300 
## 2 PR1B           1     6    0.668  0.445  
## 3 CPI8           1     6    2.55   0.161  
## 4 CAB            1     6    6.52   0.0433 
## 5 BGLU2          1     6    1.96   0.211  
## 6 HSP70          1     6    0.0221 0.887  
## 7 PTI5           1     6   24.4    0.00260
## 8 13-LOX         1     6    4.36   0.0818 
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 RBOHD          1     6    8.01   0.0300 
## 2 PR1B           1     6    0.668  0.445  
## 3 CPI8           1     6    2.55   0.161  
## 4 CAB            1     6    6.52   0.0433 
## 5 BGLU2          1     6    1.96   0.211  
## 6 HSP70          1     6    0.0221 0.887  
## 7 PTI5           1     6   24.4    0.00260
## 8 13-LOX         1     6    4.36   0.0818 
## Wilcoxon effect size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.306 RBOHD          4     4 moderate 
## 2 measurement noninoculated inoculated   0.408 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated   0.510 CPI8           4     4 large    
## 4 measurement noninoculated inoculated   0.816 CAB            4     4 large    
## 5 measurement noninoculated inoculated   0.408 BGLU2          4     4 moderate 
## 6 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.306 RBOHD          4     4 moderate 
## 2 scaled noninoculated inoculated   0.408 PR1B           4     4 moderate 
## 3 scaled noninoculated inoculated   0.510 CPI8           4     4 large    
## 4 scaled noninoculated inoculated   0.816 CAB            4     4 large    
## 5 scaled noninoculated inoculated   0.408 BGLU2          4     4 moderate 
## 6 scaled noninoculated inoculated   0.204 HSP70          4     4 small    
## 7 scaled noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 scaled noninoculated inoculated   0.816 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.928 RBOHD          4     4 large    
## 2 measurement noninoculated inoculated  -0.769 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated  -1.44  CPI8           4     4 large    
## 4 measurement noninoculated inoculated  -1.94  CAB            4     4 large    
## 5 measurement noninoculated inoculated   0.906 BGLU2          4     4 large    
## 6 measurement noninoculated inoculated  -0.560 HSP70          4     4 moderate 
## 7 measurement noninoculated inoculated  -2.12  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -2.81  13-LOX         4     4 large    
## # A tibble: 8 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated  -0.928 RBOHD          4     4 large    
## 2 scaled noninoculated inoculated  -0.769 PR1B           4     4 moderate 
## 3 scaled noninoculated inoculated  -1.44  CPI8           4     4 large    
## 4 scaled noninoculated inoculated  -1.94  CAB            4     4 large    
## 5 scaled noninoculated inoculated   0.906 BGLU2          4     4 large    
## 6 scaled noninoculated inoculated  -0.560 HSP70          4     4 moderate 
## 7 scaled noninoculated inoculated  -2.12  PTI5           4     4 large    
## 8 scaled noninoculated inoculated  -2.81  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p", colnames(result$stat.test))]
print(res)
##   transcript        group1     group2          p      p.adj p.adj.signif
## 1      RBOHD noninoculated inoculated 0.20000000 0.26666667           ns
## 2       PR1B noninoculated inoculated 0.28571429 0.32653061           ns
## 3       CPI8 noninoculated inoculated 0.12857143 0.25714286           ns
## 4        CAB noninoculated inoculated 0.01428571 0.03809524            *
## 5      BGLU2 noninoculated inoculated 0.20000000 0.26666667           ns
## 6      HSP70 noninoculated inoculated 0.60000000 0.60000000           ns
## 7       PTI5 noninoculated inoculated 0.01428571 0.03809524            *
## 8     13-LOX noninoculated inoculated 0.01428571 0.03809524            *
output_text = capture.output(paste0("Test - roots - 218"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
output_text = capture.output(res)
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")

ggsave(
  filename = file.path(my_dir, "Test_roots.218_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_roots.218_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)
keep = c("magrittr", "ggplot2", "%!in%", "pal", "my_dir", "fr")
all_objs = ls()
# Objects that are functions or named in keep
keep = c(lsf.str(), intersect(all_objs, keep))
rm(list = setdiff(all_objs, keep))

5 Exp1

5.1 input

fp = file.path('..', "input")
list.files(fp)
## [1] "README.md"     "Table S1.xlsx" "Table S2.xlsx" "Table S7.xlsx"
fn = 'Table S7.xlsx'

PS218 = openxlsx::read.xlsx(xlsxFile = file.path(fp, fn),
                            sheet = 'Exp1',
                            startRow = 1,
                            colNames = TRUE,
                            rowNames = FALSE,
                            detectDates = FALSE,
                            skipEmptyRows = TRUE,
                            skipEmptyCols = TRUE,
                            rows = NULL,
                            cols = NULL,
                            check.names = FALSE,
                            sep.names = ".",
                            namedRegion = NULL,
                            na.strings = "NA",
                            fillMergedCells = FALSE)

data.table::setDT(PS218)
PS218[, Sample.ID := NULL]
PS218[, Genotype := NULL]
PS218[, Strain := NULL]
PS218$Tissue = as.factor(trimws(PS218$Tissue))
PS218$Treatment = factor(trimws(PS218$Treatment), levels = c("noninoculated", "inoculated"))

table(PS218$Time)
## 
## 2 h 4 h 6 h 
##  16  16  16

5.2 params

groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("BGLU2", "HSP70", "PTI5", "13-LOX")

5.3 results

run_analysis_for_time <- function(data
                                  , time_point
                                  , myorder
                                  , pal
                                  , groupvars
                                  , measurevar
                                  , y_scales1
                                  , y_scales2
                                  , y_scales3
                                  , y_scales4
                                  , my_dir
                                  , plot_gene_expression
                                  , test_and_plot) {

  temp = data[data$Time == time_point, ]
  temp[, Time := NULL]

  shoots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
  roots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$roots
  

  result_shoots = test_and_plot(data_long_raw = shoots,
                                myorder = myorder,
                                pal = pal,
                                what = paste("shoots", time_point),
                                plot_gene_expression_func = plot_gene_expression,
                                groupvars = groupvars,
                                y_scales1 = y_scales1,
                                y_scales2 = y_scales2)
  res = result_shoots$stat.test[, grep("transcript|group1|group2|^p", colnames(result_shoots$stat.test))]
  print(res)
  
  output_text = capture.output(paste0("Exp 1 shoots - ", time_point))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  output_text = capture.output(res)
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Exp1_shoots.", gsub(" ", "", time_point), "_1.pdf")),
         plot = result_shoots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Exp1_shoots.", gsub(" ", "", time_point), "_2.pdf")),
         plot = result_shoots$plot2, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  

  result_roots = test_and_plot(data_long_raw = roots,
                               myorder = myorder,
                               pal = pal,
                               what = paste("roots", time_point),
                               plot_gene_expression_func = plot_gene_expression,
                               groupvars = groupvars,
                               y_scales1 = y_scales3,
                               y_scales2 = y_scales4)
  res = result_roots$stat.test[, grep("transcript|group1|group2|^p", colnames(result_roots$stat.test))]
  print(res)
  
  output_text = capture.output(paste0("Exp 1 roots - ", time_point))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  output_text = capture.output(res)
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Exp1_roots.", gsub(" ", ".", time_point), "_1.pdf")),
         plot = result_roots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Exp1_roots.", gsub(" ", ".", time_point), "_2.pdf")),
         plot = result_roots$plot2, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  

  list(shoots = result_shoots, roots = result_roots)
}


y_scales1 = list(transcript == "PTI5" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 1))
                 , transcript == "13-LOX" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 1)))
y_scales2 = list(transcript == "BGLU2" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 , transcript == "HSP70" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 )
y_scales3 = list(transcript == "PTI5" ~ scale_y_continuous(limits = c(0, 11), breaks = seq(0, 11, 1))
                 , transcript == "13-LOX" ~ scale_y_continuous(limits = c(0, 11), breaks = seq(0, 11, 1)))
y_scales4 = list(transcript == "BGLU2" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 , transcript == "HSP70" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 )

results_2h = run_analysis_for_time(data = PS218
                                   , time_point = "2 h"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , y_scales1
                                   , y_scales2
                                   , y_scales3
                                   , y_scales4
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.907   0.335     1
## 2 HSP70      0.954   0.749     1
## 3 PTI5       0.931   0.524     1
## 4 13-LOX     0.956   0.773     1
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.907   0.335     1
## 2 HSP70      0.954   0.749     1
## 3 PTI5       0.931   0.524     1
## 4 13-LOX     0.956   0.773     1
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6    1.91   0.216
## 2 HSP70          1     6    0.0663 0.805
## 3 PTI5           1     6    0.156  0.707
## 4 13-LOX         1     6    0.0281 0.872
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6    1.91   0.216
## 2 HSP70          1     6    0.0663 0.805
## 3 PTI5           1     6    0.156  0.707
## 4 13-LOX         1     6    0.0281 0.872
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.510 BGLU2          4     4 large    
## 2 measurement noninoculated inoculated   0.102 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.612 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0     13-LOX         4     4 small    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.510 BGLU2          4     4 large    
## 2 scaled noninoculated inoculated   0.102 HSP70          4     4 small    
## 3 scaled noninoculated inoculated   0.612 PTI5           4     4 large    
## 4 scaled noninoculated inoculated   0     13-LOX         4     4 small    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  1.13   BGLU2          4     4 large     
## 2 measurement noninoculated inoculated -0.176  HSP70          4     4 negligible
## 3 measurement noninoculated inoculated -1.19   PTI5           4     4 large     
## 4 measurement noninoculated inoculated  0.0435 13-LOX         4     4 negligible
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 scaled noninoculated inoculated  1.13   BGLU2          4     4 large     
## 2 scaled noninoculated inoculated -0.176  HSP70          4     4 negligible
## 3 scaled noninoculated inoculated -1.19   PTI5           4     4 large     
## 4 scaled noninoculated inoculated  0.0435 13-LOX         4     4 negligible

##   transcript        group1     group2         p     p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.1000000 0.2285714           ns
## 2      HSP70 noninoculated inoculated 0.8428571 0.9571429           ns
## 3       PTI5 noninoculated inoculated 0.1142857 0.2285714           ns
## 4     13-LOX noninoculated inoculated 0.9571429 0.9571429           ns

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.886  0.252  0.252 
## 2 HSP70      0.787  0.0209 0.0626
## 3 PTI5       0.814  0.0404 0.0808
## 4 13-LOX     0.761  0.0108 0.0434
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.886  0.252  0.252 
## 2 HSP70      0.787  0.0209 0.0626
## 3 PTI5       0.814  0.0404 0.0808
## 4 13-LOX     0.761  0.0108 0.0434
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     2.63  0.166   
## 2 HSP70          1     6    38.9   0.000785
## 3 PTI5           1     6     4.06  0.0906  
## 4 13-LOX         1     6     0.370 0.566   
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     2.63  0.166   
## 2 HSP70          1     6    38.9   0.000785
## 3 PTI5           1     6     4.06  0.0906  
## 4 13-LOX         1     6     0.370 0.566   
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.267 BGLU2          4     3 small    
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.267 BGLU2          4     3 small    
## 2 scaled noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 scaled noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 scaled noninoculated inoculated   0.816 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.789 BGLU2          4     3 moderate 
## 2 measurement noninoculated inoculated  -0.823 HSP70          4     4 large    
## 3 measurement noninoculated inoculated  -1.99  PTI5           4     4 large    
## 4 measurement noninoculated inoculated -11.7   13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated  -0.789 BGLU2          4     3 moderate 
## 2 scaled noninoculated inoculated  -0.823 HSP70          4     4 large    
## 3 scaled noninoculated inoculated  -1.99  PTI5           4     4 large    
## 4 scaled noninoculated inoculated -11.7   13-LOX         4     4 large

##   transcript        group1     group2          p      p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.40000000 0.40000000           ns
## 2      HSP70 noninoculated inoculated 0.35714286 0.40000000           ns
## 3       PTI5 noninoculated inoculated 0.01428571 0.02857143            *
## 4     13-LOX noninoculated inoculated 0.01428571 0.02857143            *

y_scales1 = list(transcript == "PTI5" ~ scale_y_continuous(limits = c(0, 5), breaks = seq(0, 11, 1))
                 , transcript == "13-LOX" ~ scale_y_continuous(limits = c(0, 5), breaks = seq(0, 11, 1)))
results_4h = run_analysis_for_time(data = PS218
                                   , time_point = "4 h"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , y_scales1
                                   , y_scales2
                                   , y_scales3
                                   , y_scales4
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.908 0.338   1     
## 2 HSP70      0.914 0.387   1     
## 3 PTI5       0.740 0.00989 0.0396
## 4 13-LOX     0.919 0.424   1     
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.908 0.338   1     
## 2 HSP70      0.914 0.387   1     
## 3 PTI5       0.740 0.00989 0.0396
## 4 13-LOX     0.919 0.424   1     
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    7.26   0.0358
## 2 HSP70          1     6    2.06   0.201 
## 3 PTI5           1     5    3.62   0.115 
## 4 13-LOX         1     6    0.0744 0.794 
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    7.26   0.0358
## 2 HSP70          1     6    2.06   0.201 
## 3 PTI5           1     5    3.62   0.115 
## 4 13-LOX         1     6    0.0744 0.794 
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.408 BGLU2          4     4 moderate 
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.668 PTI5           4     3 large    
## 4 measurement noninoculated inoculated   0.510 13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.408 BGLU2          4     4 moderate 
## 2 scaled noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 scaled noninoculated inoculated   0.668 PTI5           4     3 large    
## 4 scaled noninoculated inoculated   0.510 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  1.14   BGLU2          4     4 large     
## 2 measurement noninoculated inoculated -0.0196 HSP70          4     4 negligible
## 3 measurement noninoculated inoculated -1.41   PTI5           4     3 large     
## 4 measurement noninoculated inoculated -0.802  13-LOX         4     4 large     
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 scaled noninoculated inoculated  1.14   BGLU2          4     4 large     
## 2 scaled noninoculated inoculated -0.0196 HSP70          4     4 negligible
## 3 scaled noninoculated inoculated -1.41   PTI5           4     3 large     
## 4 scaled noninoculated inoculated -0.802  13-LOX         4     4 large

##   transcript        group1     group2          p     p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.20000000 0.4000000           ns
## 2      HSP70 noninoculated inoculated 0.97142857 0.9714286           ns
## 3       PTI5 noninoculated inoculated 0.04285714 0.1714286           ns
## 4     13-LOX noninoculated inoculated 0.37142857 0.4952381           ns

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.908  0.343  0.832
## 2 HSP70      0.904  0.315  0.832
## 3 PTI5       0.781  0.0263 0.105
## 4 13-LOX     0.898  0.277  0.832
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.908  0.343  0.832
## 2 HSP70      0.904  0.315  0.832
## 3 PTI5       0.781  0.0263 0.105
## 4 13-LOX     0.898  0.277  0.832
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    0.0398 0.849 
## 2 HSP70          1     6    0.494  0.508 
## 3 PTI5           1     5   13.7    0.0140
## 4 13-LOX         1     6    2.52   0.164 
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    0.0398 0.849 
## 2 HSP70          1     6    0.494  0.508 
## 3 PTI5           1     5   13.7    0.0140
## 4 13-LOX         1     6    2.52   0.164 
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.204 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0.714 HSP70          4     4 large    
## 3 measurement noninoculated inoculated   0.802 PTI5           3     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.204 BGLU2          4     4 small    
## 2 scaled noninoculated inoculated   0.714 HSP70          4     4 large    
## 3 scaled noninoculated inoculated   0.802 PTI5           3     4 large    
## 4 scaled noninoculated inoculated   0.816 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated   0.101 BGLU2          4     4 negligible
## 2 measurement noninoculated inoculated  -2.66  HSP70          4     4 large     
## 3 measurement noninoculated inoculated  -1.62  PTI5           3     4 large     
## 4 measurement noninoculated inoculated  -2.65  13-LOX         4     4 large     
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 scaled noninoculated inoculated   0.101 BGLU2          4     4 negligible
## 2 scaled noninoculated inoculated  -2.66  HSP70          4     4 large     
## 3 scaled noninoculated inoculated  -1.62  PTI5           3     4 large     
## 4 scaled noninoculated inoculated  -2.65  13-LOX         4     4 large

##   transcript        group1     group2          p      p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.77142857 0.77142857           ns
## 2      HSP70 noninoculated inoculated 0.01428571 0.02857143            *
## 3       PTI5 noninoculated inoculated 0.04285714 0.05714286           ns
## 4     13-LOX noninoculated inoculated 0.01428571 0.02857143            *

results_6h = run_analysis_for_time(data = PS218
                                   , time_point = "6 h"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , y_scales1
                                   , y_scales2
                                   , y_scales3
                                   , y_scales4
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.853   0.103 0.410
## 2 HSP70      0.975   0.932 1    
## 3 PTI5       0.891   0.239 0.718
## 4 13-LOX     0.958   0.793 1    
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.853   0.103 0.410
## 2 HSP70      0.975   0.932 1    
## 3 PTI5       0.891   0.239 0.718
## 4 13-LOX     0.958   0.793 1    
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6     9.33  0.0224
## 2 HSP70          1     6     1.76  0.233 
## 3 PTI5           1     6     0.632 0.457 
## 4 13-LOX         1     6     0.157 0.705 
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6     9.33  0.0224
## 2 HSP70          1     6     1.76  0.233 
## 3 PTI5           1     6     0.632 0.457 
## 4 13-LOX         1     6     0.157 0.705 
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.204 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.714 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.510 13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.204 BGLU2          4     4 small    
## 2 scaled noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 scaled noninoculated inoculated   0.714 PTI5           4     4 large    
## 4 scaled noninoculated inoculated   0.510 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  -0.783 BGLU2          4     4 moderate  
## 2 measurement noninoculated inoculated  -0.167 HSP70          4     4 negligible
## 3 measurement noninoculated inoculated  -2.26  PTI5           4     4 large     
## 4 measurement noninoculated inoculated  -1.01  13-LOX         4     4 large     
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 scaled noninoculated inoculated  -0.783 BGLU2          4     4 moderate  
## 2 scaled noninoculated inoculated  -0.167 HSP70          4     4 negligible
## 3 scaled noninoculated inoculated  -2.26  PTI5           4     4 large     
## 4 scaled noninoculated inoculated  -1.01  13-LOX         4     4 large

##   transcript        group1     group2          p     p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.41428571 0.5523810           ns
## 2      HSP70 noninoculated inoculated 0.84285714 0.8428571           ns
## 3       PTI5 noninoculated inoculated 0.04285714 0.1714286           ns
## 4     13-LOX noninoculated inoculated 0.14285714 0.2857143           ns

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.854  0.105  0.315
## 2 HSP70      0.930  0.514  0.514
## 3 PTI5       0.839  0.0740 0.296
## 4 13-LOX     0.861  0.122  0.315
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.854  0.105  0.315
## 2 HSP70      0.930  0.514  0.514
## 3 PTI5       0.839  0.0740 0.296
## 4 13-LOX     0.861  0.122  0.315
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    1.48   0.269 
## 2 HSP70          1     6    0.0375 0.853 
## 3 PTI5           1     6    1.80   0.228 
## 4 13-LOX         1     6    6.91   0.0391
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    1.48   0.269 
## 2 HSP70          1     6    0.0375 0.853 
## 3 PTI5           1     6    1.80   0.228 
## 4 13-LOX         1     6    6.91   0.0391
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.714 BGLU2          4     4 large    
## 2 measurement noninoculated inoculated   0.306 HSP70          4     4 moderate 
## 3 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.714 BGLU2          4     4 large    
## 2 scaled noninoculated inoculated   0.306 HSP70          4     4 moderate 
## 3 scaled noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 scaled noninoculated inoculated   0.816 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -2.10  BGLU2          4     4 large    
## 2 measurement noninoculated inoculated  -0.886 HSP70          4     4 large    
## 3 measurement noninoculated inoculated  -4.39  PTI5           4     4 large    
## 4 measurement noninoculated inoculated  -2.64  13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated  -2.10  BGLU2          4     4 large    
## 2 scaled noninoculated inoculated  -0.886 HSP70          4     4 large    
## 3 scaled noninoculated inoculated  -4.39  PTI5           4     4 large    
## 4 scaled noninoculated inoculated  -2.64  13-LOX         4     4 large

##   transcript        group1     group2          p      p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.01428571 0.01904762            *
## 2      HSP70 noninoculated inoculated 0.27142857 0.27142857           ns
## 3       PTI5 noninoculated inoculated 0.01428571 0.01904762            *
## 4     13-LOX noninoculated inoculated 0.01428571 0.01904762            *

keep = c("magrittr", "ggplot2", "%!in%", "pal", "my_dir", "fr")
all_objs = ls()
# Objects that are functions or named in keep
keep = c(lsf.str(), intersect(all_objs, keep))
rm(list = setdiff(all_objs, keep))

6 Exp2

6.1 input

fp = file.path('..', "input")
list.files(fp)
## [1] "README.md"     "Table S1.xlsx" "Table S2.xlsx" "Table S7.xlsx"
fn = 'Table S7.xlsx'

PS218 = openxlsx::read.xlsx(xlsxFile = file.path(fp, fn),
                            sheet = 'Exp2',
                            startRow = 1,
                            colNames = TRUE,
                            rowNames = FALSE,
                            detectDates = FALSE,
                            skipEmptyRows = TRUE,
                            skipEmptyCols = TRUE,
                            rows = NULL,
                            cols = NULL,
                            check.names = FALSE,
                            sep.names = ".",
                            namedRegion = NULL,
                            na.strings = "NA",
                            fillMergedCells = FALSE)

data.table::setDT(PS218)
PS218[, Sample.ID := NULL]
PS218[, Genotype := NULL]
PS218[, Strain := NULL]
PS218$Tissue = as.factor(trimws(PS218$Tissue))
PS218$Treatment = factor(trimws(PS218$Treatment), levels = c("noninoculated", "inoculated"))

table(PS218$Time)
## 
##    1 h  1 min 15 min    2 h 30 min 
##      8      8      8      8      8

6.2 params

groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("BGLU2", "HSP70", "PTI5", "13-LOX")

6.3 results

run_analysis_for_time <- function(data
                                  , time_point
                                  , myorder
                                  , pal
                                  , groupvars
                                  , measurevar
                                  , y_scales1
                                  , y_scales2
                                  , my_dir
                                  , plot_gene_expression
                                  , test_and_plot
                                  ) {

  temp = data[data$Time == time_point, ]
  temp[, Time := NULL]

  roots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$roots
  

  result_roots = test_and_plot(data_long_raw = roots,
                               myorder = myorder,
                               pal = pal,
                               what = paste("roots", time_point),
                               plot_gene_expression_func = plot_gene_expression,
                               groupvars = groupvars,
                               y_scales1,
                               y_scales2)
  res = result_roots$stat.test[, grep("transcript|group1|group2|^p", colnames(result_roots$stat.test))]
  print(res)
  
  output_text = capture.output(paste0("Exp 2 roots - ", time_point))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  output_text = capture.output(res)
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")


  
  ggsave(filename = file.path(my_dir, paste0("Exp2_roots.", gsub(" ", ".", time_point), "_1.pdf")),
         plot = result_roots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Exp2_roots.", gsub(" ", ".", time_point), "_2.pdf")),
         plot = result_roots$plot2, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  

  list(roots = result_roots)
}


y_scales1 = list(transcript == "PTI5" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 1))
                 , transcript == "13-LOX" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 1)))
y_scales2 = list(transcript == "BGLU2" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 , transcript == "HSP70" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 )

results_1min = run_analysis_for_time(data = PS218
                                     , time_point = "1 min"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.898   0.278 0.835
## 2 HSP70      0.956   0.774 0.863
## 3 PTI5       0.920   0.431 0.863
## 4 13-LOX     0.884   0.205 0.819
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.898   0.278 0.835
## 2 HSP70      0.956   0.774 0.863
## 3 PTI5       0.920   0.431 0.863
## 4 13-LOX     0.884   0.205 0.819
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6  1.63     0.249
## 2 HSP70          1     6  0.530    0.494
## 3 PTI5           1     6  0.000531 0.982
## 4 13-LOX         1     6  0.652    0.450
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6  1.63     0.249
## 2 HSP70          1     6  0.530    0.494
## 3 PTI5           1     6  0.000531 0.982
## 4 13-LOX         1     6  0.652    0.450
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.612 BGLU2          4     4 large    
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.102 PTI5           4     4 small    
## 4 measurement noninoculated inoculated   0.204 13-LOX         4     4 small    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.612 BGLU2          4     4 large    
## 2 scaled noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 scaled noninoculated inoculated   0.102 PTI5           4     4 small    
## 4 scaled noninoculated inoculated   0.204 13-LOX         4     4 small    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  -1.05  BGLU2          4     4 large     
## 2 measurement noninoculated inoculated  -0.798 HSP70          4     4 moderate  
## 3 measurement noninoculated inoculated   0.430 PTI5           4     4 small     
## 4 measurement noninoculated inoculated   0.137 13-LOX         4     4 negligible
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 scaled noninoculated inoculated  -1.05  BGLU2          4     4 large     
## 2 scaled noninoculated inoculated  -0.798 HSP70          4     4 moderate  
## 3 scaled noninoculated inoculated   0.430 PTI5           4     4 small     
## 4 scaled noninoculated inoculated   0.137 13-LOX         4     4 negligible

##   transcript        group1     group2         p     p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.1142857 0.4571429           ns
## 2      HSP70 noninoculated inoculated 0.2714286 0.5428571           ns
## 3       PTI5 noninoculated inoculated 0.5857143 0.7809524           ns
## 4     13-LOX noninoculated inoculated 0.8000000 0.8000000           ns

results_15min = run_analysis_for_time(data = PS218
                                     , time_point = "15 min"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.860   0.120 0.480
## 2 HSP70      0.967   0.874 1    
## 3 PTI5       0.973   0.920 1    
## 4 13-LOX     0.937   0.578 1    
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.860   0.120 0.480
## 2 HSP70      0.967   0.874 1    
## 3 PTI5       0.973   0.920 1    
## 4 13-LOX     0.937   0.578 1    
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     2.83  0.143
## 2 HSP70          1     6     0.795 0.407
## 3 PTI5           1     6     0.890 0.382
## 4 13-LOX         1     6     3.67  0.104
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     2.83  0.143
## 2 HSP70          1     6     0.795 0.407
## 3 PTI5           1     6     0.890 0.382
## 4 13-LOX         1     6     3.67  0.104
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.306 BGLU2          4     4 moderate 
## 2 measurement noninoculated inoculated   0.816 HSP70          4     4 large    
## 3 measurement noninoculated inoculated   0.408 PTI5           4     4 moderate 
## 4 measurement noninoculated inoculated   0.204 13-LOX         4     4 small    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.306 BGLU2          4     4 moderate 
## 2 scaled noninoculated inoculated   0.816 HSP70          4     4 large    
## 3 scaled noninoculated inoculated   0.408 PTI5           4     4 moderate 
## 4 scaled noninoculated inoculated   0.204 13-LOX         4     4 small    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.662 BGLU2          4     4 moderate 
## 2 measurement noninoculated inoculated   2.07  HSP70          4     4 large    
## 3 measurement noninoculated inoculated  -0.883 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.623 13-LOX         4     4 moderate 
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.662 BGLU2          4     4 moderate 
## 2 scaled noninoculated inoculated   2.07  HSP70          4     4 large    
## 3 scaled noninoculated inoculated  -0.883 PTI5           4     4 large    
## 4 scaled noninoculated inoculated   0.623 13-LOX         4     4 moderate

##   transcript        group1     group2          p      p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.30000000 0.30000000           ns
## 2      HSP70 noninoculated inoculated 0.01428571 0.05714286           ns
## 3       PTI5 noninoculated inoculated 0.22857143 0.30000000           ns
## 4     13-LOX noninoculated inoculated 0.22857143 0.30000000           ns

results_30min = run_analysis_for_time(data = PS218
                                     , time_point = "30 min"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.922   0.445     1
## 2 HSP70      0.934   0.550     1
## 3 PTI5       0.917   0.407     1
## 4 13-LOX     0.976   0.937     1
## # A tibble: 4 × 4
##   transcript     W p.value p.adj
##   <fct>      <dbl>   <dbl> <dbl>
## 1 BGLU2      0.922   0.445     1
## 2 HSP70      0.934   0.550     1
## 3 PTI5       0.917   0.407     1
## 4 13-LOX     0.976   0.937     1
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.752 0.419
## 2 HSP70          1     6     1.64  0.248
## 3 PTI5           1     6     0.310 0.598
## 4 13-LOX         1     6     1.36  0.288
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.752 0.419
## 2 HSP70          1     6     1.64  0.248
## 3 PTI5           1     6     0.310 0.598
## 4 13-LOX         1     6     1.36  0.288
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.102 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0.408 HSP70          4     4 moderate 
## 3 measurement noninoculated inoculated   0.612 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.204 13-LOX         4     4 small    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.102 BGLU2          4     4 small    
## 2 scaled noninoculated inoculated   0.408 HSP70          4     4 moderate 
## 3 scaled noninoculated inoculated   0.612 PTI5           4     4 large    
## 4 scaled noninoculated inoculated   0.204 13-LOX         4     4 small    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated -0.0596 BGLU2          4     4 negligible
## 2 measurement noninoculated inoculated -0.800  HSP70          4     4 moderate  
## 3 measurement noninoculated inoculated -1.16   PTI5           4     4 large     
## 4 measurement noninoculated inoculated -0.497  13-LOX         4     4 small     
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 scaled noninoculated inoculated -0.0596 BGLU2          4     4 negligible
## 2 scaled noninoculated inoculated -0.800  HSP70          4     4 moderate  
## 3 scaled noninoculated inoculated -1.16   PTI5           4     4 large     
## 4 scaled noninoculated inoculated -0.497  13-LOX         4     4 small

##   transcript        group1     group2          p     p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.94285714 0.9428571           ns
## 2      HSP70 noninoculated inoculated 0.27142857 0.5428571           ns
## 3       PTI5 noninoculated inoculated 0.08571429 0.3428571           ns
## 4     13-LOX noninoculated inoculated 0.51428571 0.6857143           ns

y_scales1 = list(transcript == "PTI5" ~ scale_y_continuous(limits = c(0, 8), breaks = seq(0, 11, 1))
                 , transcript == "13-LOX" ~ scale_y_continuous(limits = c(0, 8), breaks = seq(0, 11, 1)))
y_scales2 = list(transcript == "BGLU2" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 , transcript == "HSP70" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 )


results_1h = run_analysis_for_time(data = PS218
                                     , time_point = "1 h"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.894  0.298  0.597 
## 2 HSP70      0.808  0.0486 0.146 
## 3 PTI5       0.768  0.0196 0.0782
## 4 13-LOX     0.905  0.363  0.597 
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.894  0.298  0.597 
## 2 HSP70      0.808  0.0486 0.146 
## 3 PTI5       0.768  0.0196 0.0782
## 4 13-LOX     0.905  0.363  0.597 
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     93.6  0.000200
## 2 HSP70          1     5      5.10 0.0736  
## 3 PTI5           1     5      1.07 0.348   
## 4 13-LOX         1     5      1.85 0.232   
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     93.6  0.000200
## 2 HSP70          1     5      5.10 0.0736  
## 3 PTI5           1     5      1.07 0.348   
## 4 13-LOX         1     5      1.85 0.232   
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0     BGLU2          3     4 small    
## 2 measurement noninoculated inoculated   0.401 HSP70          3     4 moderate 
## 3 measurement noninoculated inoculated   0.668 PTI5           3     4 large    
## 4 measurement noninoculated inoculated   0.802 13-LOX         3     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0     BGLU2          3     4 small    
## 2 scaled noninoculated inoculated   0.401 HSP70          3     4 moderate 
## 3 scaled noninoculated inoculated   0.668 PTI5           3     4 large    
## 4 scaled noninoculated inoculated   0.802 13-LOX         3     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.393 BGLU2          3     4 small    
## 2 measurement noninoculated inoculated  -0.912 HSP70          3     4 large    
## 3 measurement noninoculated inoculated  -1.35  PTI5           3     4 large    
## 4 measurement noninoculated inoculated  -3.09  13-LOX         3     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated  -0.393 BGLU2          3     4 small    
## 2 scaled noninoculated inoculated  -0.912 HSP70          3     4 large    
## 3 scaled noninoculated inoculated  -1.35  PTI5           3     4 large    
## 4 scaled noninoculated inoculated  -3.09  13-LOX         3     4 large

##   transcript        group1     group2          p     p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.40000000 0.4000000           ns
## 2      HSP70 noninoculated inoculated 0.34285714 0.4000000           ns
## 3       PTI5 noninoculated inoculated 0.07142857 0.1428571           ns
## 4     13-LOX noninoculated inoculated 0.04285714 0.1428571           ns

y_scales1 = list(transcript == "PTI5" ~ scale_y_continuous(limits = c(0, 12), breaks = seq(0, 11, 1))
                 , transcript == "13-LOX" ~ scale_y_continuous(limits = c(0, 12), breaks = seq(0, 11, 1)))
y_scales2 = list(transcript == "BGLU2" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 , transcript == "HSP70" ~ scale_y_continuous(limits = c(0, 3), breaks = seq(0, 11, 0.5))
                 )

results_2h = run_analysis_for_time(data = PS218
                                     , time_point = "2 h"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## Shapiro-Wilk Normality Test
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.886  0.215  0.215 
## 2 HSP70      0.760  0.0106 0.0422
## 3 PTI5       0.838  0.0720 0.144 
## 4 13-LOX     0.772  0.0144 0.0431
## # A tibble: 4 × 4
##   transcript     W p.value  p.adj
##   <fct>      <dbl>   <dbl>  <dbl>
## 1 BGLU2      0.886  0.215  0.215 
## 2 HSP70      0.760  0.0106 0.0422
## 3 PTI5       0.838  0.0720 0.144 
## 4 13-LOX     0.772  0.0144 0.0431
## Quantile-Quantile plots
## Levene's test for homogeneity of variance across groups
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.151 0.711
## 2 HSP70          1     6     0.926 0.373
## 3 PTI5           1     6     2.14  0.194
## 4 13-LOX         1     6     0.634 0.456
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.151 0.711
## 2 HSP70          1     6     0.926 0.373
## 3 PTI5           1     6     2.14  0.194
## 4 13-LOX         1     6     0.634 0.456
## Wilcoxon effect size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.102 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0     HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.102 BGLU2          4     4 small    
## 2 scaled noninoculated inoculated   0     HSP70          4     4 small    
## 3 scaled noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 scaled noninoculated inoculated   0.816 13-LOX         4     4 large    
## Cohen's d Measure of Effect Size
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.413 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated  -0.442 HSP70          4     4 small    
## 3 measurement noninoculated inoculated  -2.37  PTI5           4     4 large    
## 4 measurement noninoculated inoculated -12.6   13-LOX         4     4 large    
## # A tibble: 4 × 8
##   .y.    group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>  <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 scaled noninoculated inoculated   0.413 BGLU2          4     4 small    
## 2 scaled noninoculated inoculated  -0.442 HSP70          4     4 small    
## 3 scaled noninoculated inoculated  -2.37  PTI5           4     4 large    
## 4 scaled noninoculated inoculated -12.6   13-LOX         4     4 large

##   transcript        group1     group2          p      p.adj p.adj.signif
## 1      BGLU2 noninoculated inoculated 0.67142857 0.81428571           ns
## 2      HSP70 noninoculated inoculated 0.81428571 0.81428571           ns
## 3       PTI5 noninoculated inoculated 0.01428571 0.02857143            *
## 4     13-LOX noninoculated inoculated 0.01428571 0.02857143            *

sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Ljubljana
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.5.2  magrittr_2.0.3
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1      libcoin_1.0-10        dplyr_1.1.4          
##  [4] farver_2.1.2          fastmap_1.2.0         TH.data_1.1-3        
##  [7] MKdescr_0.8           rpart_4.1.24          digest_0.6.37        
## [10] lifecycle_1.0.4       arrangements_1.1.9    survival_3.8-3       
## [13] compiler_4.4.1        rlang_1.1.5           sass_0.4.10          
## [16] tools_4.4.1           utf8_1.2.5            yaml_2.3.10          
## [19] data.table_1.17.0     knitr_1.50            ggsignif_0.6.4       
## [22] labeling_0.4.3        RColorBrewer_1.1-3    multcomp_1.4-28      
## [25] abind_1.4-8           MKinfer_1.2           withr_3.0.2          
## [28] purrr_1.0.4           ggh4x_0.3.0           nnet_7.3-20          
## [31] grid_4.4.1            stats4_4.4.1          ggpubr_0.6.0         
## [34] jomo_2.7-6            mice_3.17.0           scales_1.4.0         
## [37] iterators_1.0.14      MASS_7.3-64           dichromat_2.0-0.1    
## [40] cli_3.6.3             mvtnorm_1.3-3         rmarkdown_2.29       
## [43] crayon_1.5.3          reformulas_0.4.1      generics_0.1.4       
## [46] exactRankTests_0.8-35 rstudioapi_0.17.1     DBI_1.2.3            
## [49] minqa_1.2.8           cachem_1.1.0          modeltools_0.2-24    
## [52] splines_4.4.1         parallel_4.4.1        mitools_2.4          
## [55] matrixStats_1.5.0     vctrs_0.6.5           boot_1.3-31          
## [58] glmnet_4.1-8          Matrix_1.7-1          sandwich_3.1-1       
## [61] jsonlite_2.0.0        carData_3.0-5         car_3.1-3            
## [64] mitml_0.4-5           rstatix_0.7.2         Formula_1.2-5        
## [67] foreach_1.5.2         tidyr_1.3.1           jquerylib_0.1.4      
## [70] glue_1.8.0            pan_1.9               nloptr_2.2.1         
## [73] codetools_0.2-20      cowplot_1.1.3         stringi_1.8.7        
## [76] gtable_0.3.6          shape_1.4.6.1         lme4_1.1-37          
## [79] gmp_0.7-5             tibble_3.2.1          pillar_1.10.2        
## [82] htmltools_0.5.8.1     miceadds_3.17-44      R6_2.6.1             
## [85] Rdpack_2.6.4          evaluate_1.0.3        lattice_0.22-6       
## [88] rbibutils_2.3         backports_1.5.0       openxlsx_4.2.8       
## [91] broom_1.0.8           bslib_0.9.0           Rcpp_1.0.14          
## [94] zip_2.3.2             nlme_3.1-166          coin_1.4-3           
## [97] xfun_0.52             zoo_1.8-14            pkgconfig_2.0.3